## REPLICATION 
## Business Owners and Executives as Politicians: The Effect on Public Policy
## produces & saves data frames of RD results 
## includes multilple bandwidths and alternative operationalizations of DVs (result located in Supplemental Information)



rm(list=ls(all=TRUE))



library(rdrobust)
library(rdd)
library(tidyverse)
library(lmtest)
library(sandwich)


## set working directory
setwd("")

source("replication_mayors_source.R")


# mayoral & fiscal data 
load("replication_data_business_mayors_RD.RData")

## subset to elections with one business executive
mayors_subs <- subset(mayors, !is.na(exec_margin))


## List of DVs
varlist <- c("ttl_revs", "ttl_own_revs", 
             "ttl_taxes", "ttl_sales_taxes", "ttl_prop_taxes",
             "ttl_chgs_miscrev", 
             "ttl_debt", "ttl_debt_issued", "st_debt",
             "ttl_exps", 
             "police", "fire", "admin", "waste", 
             "roads", "parks", "libraries",
             "health", "housing","welfare",
             "basic_services", "developmental", "redistributive")




#################################################
####                                         ####
#### RD ESTIMATES ACROSS MULTIPLE BANDWIDTHS ####
####                                         ####
#################################################


h <- seq(.01, .25, .01)

results <- data.frame()

for(n in 1:length(h)){
    
    mayors_subs$wgts <- (1 - abs(mayors_subs$exec_margin / h[n])) * 
        (abs(mayors_subs$exec_margin) <= h[n])
 
    
#### base
models_base <- lapply(varlist, function(x) {
    
    lm(substitute(i ~ (exec_margin > 0) * exec_margin, list(i = as.name(paste0("pc_", x, "_lead2")))), 
       data = mayors_subs, weights = wgts, subset = wgts > 0)
    
})

base_conv <- lapply(models_base, summary)

base_robust <- lapply(models_base, function(x) {
    
    coeftest(x, vcov=vcovHC(x, type="HC1")) 
})



for(x in 1:length(varlist)){
temp_results <- data.frame(cbind(as.character(varlist[[x]]), 
                                 format_num(base_conv[[x]]$coefficients[2,1]),
                                 format_num(base_robust[[x]][2,2]),
                                 format_num(base_robust[[x]][2,4]),
                                 format_num(h[n], digits = 2),
                                 format_num(base_conv[[x]]$df[2] + base_conv[[x]]$df[3], digits = 0),
                                 "base"))

names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")

results <- rbind(results, temp_results)

}}



for(n in 1:length(h)){
    
    mayors_subs$wgts <- (1 - abs(mayors_subs$exec_margin / h[n])) * 
        (abs(mayors_subs$exec_margin) <= h[n])

#### covariates
models_covs <- lapply(varlist, function(x) {
    
    lm(substitute(i ~ (exec_margin > 0) * exec_margin + 
                      j +
                      log(ipopulation) +
                      ipctwhite +
                      log(imedhhinc) +
                      log(imedhouseval),
                  list(i = as.name(paste0("pc_", x, "_lead2")), 
                       j = as.name(paste0("pc_", x, "_lag1")))), 
       data = mayors_subs, weights = wgts, subset = wgts > 0)
    
})

covs_conv <- lapply(models_covs, summary)

covs_robust <- lapply(models_covs, function(x) {
    
    coeftest(x, vcov=vcovHC(x, type="HC1")) 
})




for(x in 1:length(varlist)){
    temp_results <- data.frame(cbind(as.character(varlist[[x]]), 
                                     format_num(covs_conv[[x]]$coefficients[2,1]),
                                     format_num(covs_robust[[x]][2,2]),
                                     format_num(covs_robust[[x]][2,4]),
                                     format_num(h[n], digits = 2),
                                     format_num(covs_conv[[x]]$df[2] + covs_conv[[x]]$df[3], digits = 0),
                                     "covs"))

    names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
    results <- rbind(results, temp_results)
}}



for(n in 1:length(h)){
    
    mayors_subs$wgts <- (1 - abs(mayors_subs$exec_margin / h[n])) * 
        (abs(mayors_subs$exec_margin) <= h[n])

#### log
models_log <- lapply(varlist, function(x) {
    
    lm(substitute(log(i + 1) ~ (exec_margin > 0) * exec_margin, list(i = as.name(paste0("pc_", x, "_lead2")))), 
       data = mayors_subs, weights = wgts, subset = wgts > 0)
    
})

log_conv <- lapply(models_log, summary)

log_robust <- lapply(models_log, function(x) {
    
    coeftest(x, vcov=vcovHC(x, type="HC1")) 
})




for(x in 1:length(varlist)){
    temp_results <- data.frame(cbind(as.character(varlist[[x]]), 
                                     format_num(log_conv[[x]]$coefficients[2,1]),
                                     format_num(log_robust[[x]][2,2]),
                                     format_num(log_robust[[x]][2,4]),
                                     format_num(h[n], digits = 2),
                                     format_num(log_conv[[x]]$df[2] + log_conv[[x]]$df[3], digits = 0),
                                     "log"))

    names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
    
    results <- rbind(results, temp_results)
}}



for(n in 1:length(h)){
    
    mayors_subs$wgts <- (1 - abs(mayors_subs$exec_margin / h[n])) * 
        (abs(mayors_subs$exec_margin) <= h[n])


#### log w/ covariates
models_log_covs <- lapply(varlist, function(x) {
    
    lm(substitute(log(i + 1) ~ (exec_margin > 0) * exec_margin + 
                      log(j + 1) +
                      log(ipopulation) +
                      ipctwhite +
                      log(imedhhinc) +
                      log(imedhouseval),
                  list(i = as.name(paste0("pc_", x, "_lead2")), 
                       j = as.name(paste0("pc_", x, "_lag1")))), 
       data = mayors_subs, weights = wgts, subset = wgts > 0)
    
})

log_covs_conv <- lapply(models_log_covs, summary)

log_covs_robust <- lapply(models_log_covs, function(x) {
    
    coeftest(x, vcov=vcovHC(x, type="HC1")) 
})



for(x in 1:length(varlist)){
    temp_results <- data.frame(cbind(as.character(varlist[[x]]), 
                                     format_num(log_covs_conv[[x]]$coefficients[2,1]),
                                     format_num(log_covs_robust[[x]][2,2]),
                                     format_num(log_covs_robust[[x]][2,4]),
                                     format_num(h[n], digits = 2),
                                     format_num(log_covs_conv[[x]]$df[2] + log_covs_conv[[x]]$df[3], digits = 0),
                                     "log covs"))

    names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
    
    results <- rbind(results, temp_results)
}}


for(n in 1:length(h)){
    
    mayors_subs$wgts <- (1 - abs(mayors_subs$exec_margin / h[n])) * 
        (abs(mayors_subs$exec_margin) <= h[n])
    

#### differenced 
models_diff <- lapply(varlist, function(x) {
    
    lm(substitute((i - j) ~ (exec_margin > 0) * exec_margin,
                  list(i = as.name(paste0("pc_", x, "_lead2")), 
                       j = as.name(paste0("pc_", x, "_lag1")))), 
       data = mayors_subs, weights = wgts, subset = wgts > 0)
    
})


diff_conv <- lapply(models_diff, summary)

diff_robust <- lapply(models_diff, function(x) {
    
    coeftest(x, vcov=vcovHC(x, type="HC1")) 
})




for(x in 1:length(varlist)){
    temp_results <- data.frame(cbind(as.character(varlist[[x]]), 
                                     format_num(diff_conv[[x]]$coefficients[2,1]),
                                     format_num(diff_robust[[x]][2,2]),
                                     format_num(diff_robust[[x]][2,4]),
                                     format_num(h[n], digits = 2),
                                     format_num(diff_conv[[x]]$df[2] + diff_conv[[x]]$df[3], digits = 0),
                                     "diff"))

    names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
    
    results <- rbind(results, temp_results)
}}


for(n in 1:length(h)){
  
  mayors_subs$wgts <- (1 - abs(mayors_subs$exec_margin / h[n])) * 
    (abs(mayors_subs$exec_margin) <= h[n])
  
  
  #### differenced w/ covs
  models_diff_covs <- lapply(varlist, function(x) {
    
    lm(substitute((i - j) ~ (exec_margin > 0) * exec_margin +
                  log(ipopulation) +
                    ipctwhite +
                    log(imedhhinc) +
                    log(imedhouseval),
                  list(i = as.name(paste0("pc_", x, "_lead2")), 
                       j = as.name(paste0("pc_", x, "_lag1")))), 
       data = mayors_subs, weights = wgts, subset = wgts > 0)
    
  })
  
  
  diff_covs_conv <- lapply(models_diff_covs, summary)
  
  diff_covs_robust <- lapply(models_diff_covs, function(x) {
    
    coeftest(x, vcov=vcovHC(x, type="HC1")) 
  })
  

  
  for(x in 1:length(varlist)){
    temp_results <- data.frame(cbind(as.character(varlist[[x]]), 
                                     format_num(diff_covs_conv[[x]]$coefficients[2,1]),
                                     format_num(diff_covs_robust[[x]][2,2]),
                                     format_num(diff_covs_robust[[x]][2,4]),
                                     format_num(h[n], digits = 2),
                                     format_num(diff_covs_conv[[x]]$df[2] + diff_covs_conv[[x]]$df[3], digits = 0),
                                     "diff covs"))

    names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
    
    results <- rbind(results, temp_results)
  }}


results <- within(results, {
    coef <- as.numeric(as.character(coef))
    robust_se <- as.numeric(as.character(robust_se))
    h <- as.numeric(as.character(h))
})


results <- arrange(results, h, outcome, specs)
results_dollars <- results

save(results_dollars, file = "replication_multi_bandwidth_results_constant_dollars.RData")




#################################################
####                                         ####
####  RD ESTIMATES - CCT OPTIMAL BANDWIDTHS  ####
####                                         ####
#################################################

results_CCT <- data.frame()


## base
rdrobust_specs <- lapply(varlist, function(x) {
  
  substitute(rdrobust(mayors_subs$i, mayors_subs$exec_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))
  
})

rdrobust_results <- lapply(rdrobust_specs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## base model    
  h_CCT <- as.numeric(rdrobust_results[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  base <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")))), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_base <- summary(base)
  
  base_robust <- coeftest(base, vcov=vcovHC(base, type="HC1")) 
  
  
  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_base$coefficients[2,1]),
                                   format_num(base_robust[2,2]),
                                   format_num(base_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_base$df[2] + summ_base$df[3], digits = 0),
                                   "CCT base"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



## w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i, exec_margin,
                                        covs = cbind(j, 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin +
                             j +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



## logged DV
## base
rdrobust_specs <- lapply(varlist, function(x) {
  
  substitute(rdrobust(log(mayors_subs$i + 1), mayors_subs$exec_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))
  
})

rdrobust_results <- lapply(rdrobust_specs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## base model    
  h_CCT <- as.numeric(rdrobust_results[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  base <-    lm(substitute(log(i + 1) ~ (exec_margin > 0) * exec_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")))), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_base <- summary(base)
  
  base_robust <- coeftest(base, vcov=vcovHC(base, type="HC1")) 
  
 
  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_base$coefficients[2,1]),
                                   format_num(base_robust[2,2]),
                                   format_num(base_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_base$df[2] + summ_base$df[3], digits = 0),
                                   "CCT log"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




## log w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(log(i + 1), exec_margin,
                                        covs = cbind(log(j + 1), 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  covs <-    lm(substitute(log(i + 1) ~ (exec_margin > 0) * exec_margin +
                             log(j + 1) +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT log covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




## diff
rdrobust_specs_diff <- lapply(varlist, function(x) {
  
  substitute(rdrobust((mayors_subs$i - mayors_subs$j), mayors_subs$exec_margin), list(i = as.name(paste0("pc_", x, "_lead2")),
                                                                                      j = as.name(paste0("pc_", x, "_lag1"))))
  
})

rdrobust_results_diff <- lapply(rdrobust_specs_diff, eval)




## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## diff model    
  h_CCT <- as.numeric(rdrobust_results_diff[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  models_diff <-    lm(substitute((i - j) ~ (exec_margin > 0) * exec_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")),
                                                                                  j = as.name(paste0("pc_", varlist[j], "_lag1")))), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_diff <- summary(models_diff)
  
  diff_robust <- coeftest(models_diff, vcov=vcovHC(models_diff, type="HC1")) 
  
 
  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_diff$coefficients[2,1]),
                                   format_num(diff_robust[2,2]),
                                   format_num(diff_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff$df[2] + summ_diff$df[3], digits = 0),
                                   "CCT diff"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}


## diff w/ covs

rdrobust_specs_diff_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i - j, exec_margin,
                                        covs = cbind(log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_diff_covs <- lapply(rdrobust_specs_diff_covs, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_diff_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  diff_covs <-    lm(substitute((i - j) ~ (exec_margin > 0) * exec_margin +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_diff_covs <- summary(diff_covs)
  
  diff_covs_robust <- coeftest(diff_covs, vcov=vcovHC(diff_covs, type="HC1")) 
  

  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_diff_covs$coefficients[2,1]),
                                   format_num(diff_covs_robust[2,2]),
                                   format_num(diff_covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff_covs$df[2] + summ_diff_covs$df[3], digits = 0),
                                   "CCT diff covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



save(results_CCT, file = "replication_CCT_bandwidth_results_constant_dollars.RData")



#################################################
####                                         ####
####  RD ESTIMATES - CCT OPTIMAL BANDWIDTHS  ####
####          SPENDING SHARES                ####
#################################################

shares_varlist <- c("police", "fire", "admin", "waste",
                    "roads", "parks", "libraries",
                    "health", "housing","welfare")


############ CCT bandwidths 
results_spending_shares_CCT <- data.frame()


## base
rdrobust_specs <- lapply(shares_varlist, function(x) {
  
  substitute(rdrobust(mayors_subs$i, mayors_subs$exec_margin), list(i = as.name(paste0("share_", x, "_lead2"))))
  
})

rdrobust_results <- lapply(rdrobust_specs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(shares_varlist)){
  
  ## base model    
  h_CCT <- as.numeric(rdrobust_results[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  base <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin, list(i = as.name(paste0("share_", shares_varlist[j], "_lead2")))), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_base <- summary(base)
  
  base_robust <- coeftest(base, vcov=vcovHC(base, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(shares_varlist[j]), 
                                   format_num(summ_base$coefficients[2,1]),
                                   format_num(base_robust[2,2]),
                                   format_num(base_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_base$df[2] + summ_base$df[3], digits = 0),
                                   "CCT base"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_spending_shares_CCT <- rbind(results_spending_shares_CCT, temp_results)
  
}



## w/ covs

rdrobust_specs_covs <- lapply(shares_varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i, exec_margin,
                                        covs = cbind(j, 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("share_", x, "_lead2")), 
                  j = as.name(paste0("share_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(shares_varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin +
                             j +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("share_", shares_varlist[j], "_lead2")), 
                                j = as.name(paste0("share_", shares_varlist[j], "_lag1"))
                           )), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(shares_varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_spending_shares_CCT <- rbind(results_spending_shares_CCT, temp_results)
  
}





## diff
rdrobust_specs_diff <- lapply(shares_varlist, function(x) {
  
  substitute(rdrobust((mayors_subs$i - mayors_subs$j), mayors_subs$exec_margin), list(i = as.name(paste0("share_", x, "_lead2")),
                                                                                      j = as.name(paste0("share_", x, "_lag1"))))
  
})

rdrobust_results_diff <- lapply(rdrobust_specs_diff, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(shares_varlist)){
  
  ## diff model    
  h_CCT <- as.numeric(rdrobust_results_diff[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  models_diff <-    lm(substitute((i - j) ~ (exec_margin > 0) * exec_margin, list(i = as.name(paste0("share_", shares_varlist[j], "_lead2")),
                                                                                  j = as.name(paste0("share_", shares_varlist[j], "_lag1")))), 
                       data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_diff <- summary(models_diff)
  
  diff_robust <- coeftest(models_diff, vcov=vcovHC(models_diff, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(shares_varlist[j]), 
                                   format_num(summ_diff$coefficients[2,1]),
                                   format_num(diff_robust[2,2]),
                                   format_num(diff_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff$df[2] + summ_diff$df[3], digits = 0),
                                   "CCT diff"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_spending_shares_CCT <- rbind(results_spending_shares_CCT, temp_results)
  
}


## diff w/ covs

rdrobust_specs_diff_covs <- lapply(shares_varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i - j, exec_margin,
                                        covs = cbind(log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("share_", x, "_lead2")), 
                  j = as.name(paste0("share_", x, "_lag1"))
             ))
})



rdrobust_results_diff_covs <- lapply(rdrobust_specs_diff_covs, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(shares_varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_diff_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  diff_covs <-    lm(substitute((i - j) ~ (exec_margin > 0) * exec_margin +
                                  log(ipopulation) +
                                  ipctwhite +
                                  log(imedhhinc) +
                                  log(imedhouseval),
                                list(i = as.name(paste0("share_", shares_varlist[j], "_lead2")), 
                                     j = as.name(paste0("share_", shares_varlist[j], "_lag1"))
                                )), 
                     data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_diff_covs <- summary(diff_covs)
  
  diff_covs_robust <- coeftest(diff_covs, vcov=vcovHC(diff_covs, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(shares_varlist[j]), 
                                   format_num(summ_diff_covs$coefficients[2,1]),
                                   format_num(diff_covs_robust[2,2]),
                                   format_num(diff_covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff_covs$df[2] + summ_diff_covs$df[3], digits = 0),
                                   "CCT diff covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_spending_shares_CCT <- rbind(results_spending_shares_CCT, temp_results)
  
}


save(results_spending_shares_CCT, file = "replication_CCT_bandwidth_results_spending_shares.RData")


#################################################
####                                         ####
####   RD ESTIMATES ACROSS MULTIPLE YEARS    ####
####                                         ####
#################################################


lead_list <- c("_lead1", "_lead2", "_lead3", "_lead4")

results_CCT_multi_year <- data.frame()


## w/ covs
for(k in 1:length(lead_list)){
  rdrobust_specs_covs <- lapply(varlist, function(x) {
    
    substitute(with(mayors_subs, rdrobust(i, exec_margin,
                                          covs = cbind(j, 
                                                       log(ipopulation),
                                                       ipctwhite,
                                                       log(imedhhinc),
                                                       log(imedhouseval)))),
               list(i = as.name(paste0("pc_", x, lead_list[k])), 
                    j = as.name(paste0("pc_", x, "_lag1"))
               ))
  })
  
  
  
  rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)
  
  
  
  ## these results use CCT optimal bandwidth but regular SEs
  for(j in 1:length(varlist)){
    
    ## covs model    
    h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
    
    mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
    
    
    covs <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin +
                               j +
                               log(ipopulation) +
                               ipctwhite +
                               log(imedhhinc) +
                               log(imedhouseval),
                             list(i = as.name(paste0("pc_", varlist[j], lead_list[k])), 
                                  j = as.name(paste0("pc_", varlist[j], "_lag1"))
                             )), 
                  data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
    
    summ_covs <- summary(covs)
    
    covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 

    
    temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                     format_num(summ_covs$coefficients[2,1]),
                                     format_num(covs_robust[2,2]),
                                     format_num(covs_robust[2,4]),
                                     format_num(h_CCT),
                                     format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                     "CCT covs",
                                     lead_list[k]))

    names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs", "lead_years")
    
    results_CCT_multi_year <- rbind(results_CCT_multi_year, temp_results)
    
  }
}



save(results_CCT_multi_year, file = "replication_CCT_bandwidth_multi_year_results_constant_dollars.RData")



##############################
######  PARTY ANALYSES  ######
##############################

mayors_subs_party <- 
  mayors %>% 
  mutate(mayor_major_party = ifelse(mayor_party == "D" | mayor_party == "R", mayor_party, NA),
         runnerup_major_party = ifelse(runnerup_party == "D" | runnerup_party == "R", runnerup_party, NA),
         two_party_race = !is.na(mayor_major_party) & !is.na(runnerup_major_party) & mayor_major_party != runnerup_major_party,
         rep_margin = case_when(two_party_race == TRUE & mayor_major_party == "R" ~ (mayor_votes/(mayor_votes + runnerup_votes)) - .50,
                                two_party_race == TRUE & runnerup_major_party == "R" ~ (runnerup_votes/(mayor_votes + runnerup_votes)) - .50),
         one_rep_race = (mayor_party != "" & mayor_party == "R" & runnerup_party != "R") | (runnerup_party != "" & runnerup_party == "R" & mayor_party != "R"),
         one_rep_margin = case_when(one_rep_race == TRUE & business_exec_candidates == 1 & 
                                      mayor_major_party == "R" ~ (mayor_votes/(mayor_votes + runnerup_votes)) - .50,
                                    one_rep_race == TRUE & business_exec_candidates == 1 & 
                                      runnerup_major_party == "R" ~ (runnerup_votes/(mayor_votes + runnerup_votes)) - .50) )

## RD density test
DCdensity(runvar = mayors_subs_party$rep_margin
          , cutpoint = 0
          , bin = NULL
          , bw = NULL
          , verbose = TRUE
          , plot = TRUE
          , ext.out = FALSE )

## results
# Using calculated bin size:  0.009 
# Using calculated bandwidth:  0.069 
# Log difference in heights is  -0.287  with SE  0.214 
# this gives a z-stat of  -1.343 
# and a p value of  0.179 
# [1] 0.1792771


DCdensity(runvar = mayors_subs_party$one_rep_margin
          , cutpoint = 0
          , bin = NULL
          , bw = NULL
          , verbose = TRUE
          , plot = TRUE
          , ext.out = FALSE )

## results
# Using calculated bin size:  0.013 
# Using calculated bandwidth:  0.069 
# Log difference in heights is  -0.033  with SE  0.321 
# this gives a z-stat of  -0.103 
# and a p value of  0.918 
# [1] 0.9181555


###########################################
#### RD Republican - Dem vs Rep subset ####
###########################################

results_CCT <- data.frame()


## base
rdrobust_specs <- lapply(varlist, function(x) {
  
  substitute(rdrobust(mayors_subs_party$i, mayors_subs_party$rep_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))
  
})

rdrobust_results <- lapply(rdrobust_specs, eval)




## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## base model    
  h_CCT <- as.numeric(rdrobust_results[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$rep_margin, h_CCT)
  
  
  base <-    lm(substitute(i ~ (rep_margin > 0) * rep_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")))), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_base <- summary(base)
  
  base_robust <- coeftest(base, vcov=vcovHC(base, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_base$coefficients[2,1]),
                                   format_num(base_robust[2,2]),
                                   format_num(base_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_base$df[2] + summ_base$df[3], digits = 0),
                                   "CCT base"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



## w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs_party, rdrobust(i, rep_margin,
                                        covs = cbind(j, 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$rep_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (rep_margin > 0) * rep_margin +
                             j +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



## logged DV
## base
rdrobust_specs <- lapply(varlist, function(x) {
  
  substitute(rdrobust(log(mayors_subs_party$i + 1), mayors_subs_party$rep_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))
  
})

rdrobust_results <- lapply(rdrobust_specs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## base model    
  h_CCT <- as.numeric(rdrobust_results[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$rep_margin, h_CCT)
  
  
  base <-    lm(substitute(log(i + 1) ~ (rep_margin > 0) * rep_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")))), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_base <- summary(base)
  
  base_robust <- coeftest(base, vcov=vcovHC(base, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_base$coefficients[2,1]),
                                   format_num(base_robust[2,2]),
                                   format_num(base_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_base$df[2] + summ_base$df[3], digits = 0),
                                   "CCT log"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




## log w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs_party, rdrobust(log(i + 1), rep_margin,
                                        covs = cbind(log(j + 1), 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$rep_margin, h_CCT)
  
  
  covs <-    lm(substitute(log(i + 1) ~ (rep_margin > 0) * rep_margin +
                             log(j + 1) +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT log covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




## diff
rdrobust_specs_diff <- lapply(varlist, function(x) {
  
  substitute(rdrobust((mayors_subs_party$i - mayors_subs_party$j), mayors_subs_party$rep_margin), list(i = as.name(paste0("pc_", x, "_lead2")),
                                                                                     j = as.name(paste0("pc_", x, "_lag1"))))
  
})

rdrobust_results_diff <- lapply(rdrobust_specs_diff, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## diff model    
  h_CCT <- as.numeric(rdrobust_results_diff[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$rep_margin, h_CCT)
  
  
  models_diff <-    lm(substitute((i - j) ~ (rep_margin > 0) * rep_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")),
                                                                                j = as.name(paste0("pc_", varlist[j], "_lag1")))), 
                       data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_diff <- summary(models_diff)
  
  diff_robust <- coeftest(models_diff, vcov=vcovHC(models_diff, type="HC1")) 
  

  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_diff$coefficients[2,1]),
                                   format_num(diff_robust[2,2]),
                                   format_num(diff_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff$df[2] + summ_diff$df[3], digits = 0),
                                   "CCT diff"))
  

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}


## diff w/ covs

rdrobust_specs_diff_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs_party, rdrobust(i - j, rep_margin,
                                        covs = cbind(log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_diff_covs <- lapply(rdrobust_specs_diff_covs, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_diff_covs[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$rep_margin, h_CCT)
  
  
  diff_covs <-    lm(substitute((i - j) ~ (rep_margin > 0) * rep_margin +
                                  log(ipopulation) +
                                  ipctwhite +
                                  log(imedhhinc) +
                                  log(imedhouseval),
                                list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                     j = as.name(paste0("pc_", varlist[j], "_lag1"))
                                )), 
                     data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_diff_covs <- summary(diff_covs)
  
  diff_covs_robust <- coeftest(diff_covs, vcov=vcovHC(diff_covs, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_diff_covs$coefficients[2,1]),
                                   format_num(diff_covs_robust[2,2]),
                                   format_num(diff_covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff_covs$df[2] + summ_diff_covs$df[3], digits = 0),
                                   "CCT diff covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




two_party_results_CCT <- 
  results_CCT %>% 
  mutate(candidate_type = "partisan_rep")




############################################
####    RD Republican - 1 rep subset    ####
############################################

results_CCT <- data.frame()


## base
rdrobust_specs <- lapply(varlist, function(x) {
  
  substitute(rdrobust(mayors_subs_party$i, mayors_subs_party$one_rep_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))
  
})

rdrobust_results <- lapply(rdrobust_specs, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## base model    
  h_CCT <- as.numeric(rdrobust_results[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$one_rep_margin, h_CCT)
  
  
  base <-    lm(substitute(i ~ (one_rep_margin > 0) * one_rep_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")))), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_base <- summary(base)
  
  base_robust <- coeftest(base, vcov=vcovHC(base, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_base$coefficients[2,1]),
                                   format_num(base_robust[2,2]),
                                   format_num(base_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_base$df[2] + summ_base$df[3], digits = 0),
                                   "CCT base"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



## w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs_party, rdrobust(i, one_rep_margin,
                                        covs = cbind(j, 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$one_rep_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (one_rep_margin > 0) * one_rep_margin +
                             j +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



## logged DV
## base
rdrobust_specs <- lapply(varlist, function(x) {
  
  substitute(rdrobust(log(mayors_subs_party$i + 1), mayors_subs_party$one_rep_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))
  
})

rdrobust_results <- lapply(rdrobust_specs, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## base model    
  h_CCT <- as.numeric(rdrobust_results[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$one_rep_margin, h_CCT)
  
  
  base <-    lm(substitute(log(i + 1) ~ (one_rep_margin > 0) * one_rep_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")))), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_base <- summary(base)
  
  base_robust <- coeftest(base, vcov=vcovHC(base, type="HC1")) 
  

  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_base$coefficients[2,1]),
                                   format_num(base_robust[2,2]),
                                   format_num(base_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_base$df[2] + summ_base$df[3], digits = 0),
                                   "CCT log"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




## log w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs_party, rdrobust(log(i + 1), one_rep_margin,
                                        covs = cbind(log(j + 1), 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$one_rep_margin, h_CCT)
  
  
  covs <-    lm(substitute(log(i + 1) ~ (one_rep_margin > 0) * one_rep_margin +
                             log(j + 1) +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  
  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT log covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




## diff
rdrobust_specs_diff <- lapply(varlist, function(x) {
  
  substitute(rdrobust((mayors_subs_party$i - mayors_subs_party$j), mayors_subs_party$one_rep_margin), list(i = as.name(paste0("pc_", x, "_lead2")),
                                                                                         j = as.name(paste0("pc_", x, "_lag1"))))
  
})

rdrobust_results_diff <- lapply(rdrobust_specs_diff, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## diff model    
  h_CCT <- as.numeric(rdrobust_results_diff[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$one_rep_margin, h_CCT)
  
  
  models_diff <-    lm(substitute((i - j) ~ (one_rep_margin > 0) * one_rep_margin, list(i = as.name(paste0("pc_", varlist[j], "_lead2")),
                                                                                        j = as.name(paste0("pc_", varlist[j], "_lag1")))), 
                       data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_diff <- summary(models_diff)
  
  diff_robust <- coeftest(models_diff, vcov=vcovHC(models_diff, type="HC1")) 
  

  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_diff$coefficients[2,1]),
                                   format_num(diff_robust[2,2]),
                                   format_num(diff_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff$df[2] + summ_diff$df[3], digits = 0),
                                   "CCT diff"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}


## diff w/ covs

rdrobust_specs_diff_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs_party, rdrobust(i - j, one_rep_margin,
                                        covs = cbind(log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_diff_covs <- lapply(rdrobust_specs_diff_covs, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_diff_covs[[j]]$bws[1, 2])
  
  mayors_subs_party$wgts_CCT <- weights_CCT(mayors_subs_party$one_rep_margin, h_CCT)
  
  
  diff_covs <-    lm(substitute((i - j) ~ (one_rep_margin > 0) * one_rep_margin +
                                  log(ipopulation) +
                                  ipctwhite +
                                  log(imedhhinc) +
                                  log(imedhouseval),
                                list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                     j = as.name(paste0("pc_", varlist[j], "_lag1"))
                                )), 
                     data = mayors_subs_party, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_diff_covs <- summary(diff_covs)
  
  diff_covs_robust <- coeftest(diff_covs, vcov=vcovHC(diff_covs, type="HC1")) 
  
  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_diff_covs$coefficients[2,1]),
                                   format_num(diff_covs_robust[2,2]),
                                   format_num(diff_covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_diff_covs$df[2] + summ_diff_covs$df[3], digits = 0),
                                   "CCT diff covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




one_rep_results_CCT <- 
  results_CCT %>% 
  mutate(candidate_type = "one_rep")


party_results_CCT <- bind_rows(two_party_results_CCT, one_rep_results_CCT)

save(party_results_CCT, file = "replication_CCT_party_results_constant_dollars.RData")




